#----Basics----

rm(list = ls())
# Set wokring directory
setwd("~/Dropbox/Whitten_Kagalwala/TS/JOP/Acceptance/KW_Replication/ADL_DGP")


# Load required libraries
library(tidyverse)
library(haven)
library(stringr)

# Load data
data <- read_dta('adl.dta')
data$time <- c(rep(1, 25), rep(2, 25), rep(3, 25), rep(4, 25), rep(5, 25), rep(6, 25), rep(7, 25))

# Reshape data
phi1_bias <- data %>%
  select(c(time, phi_y, phi_x, matches("^bias_(.*)_phi$"))) %>%
  pivot_longer(bias_adl_phi:bias_adl33_phi, names_to = "model", values_to = "phi1_bias") %>%
  mutate(model = gsub("^bias_(.*)_phi", "\\1", model))
phi1_type1 <- data %>%
  select(c(time, phi_y, phi_x, matches("^phi_(.*)_t1$"), matches("^phi1_(.*)_t1$"))) %>%
  pivot_longer(phi_adl_t1:phi1_adl33_t1, names_to = "model", values_to = "phi1_type1") %>%
  mutate(model = gsub("^ph.*_(.*)_t1", "\\1", model))
phi1_power <- data %>%
  select(c(time, phi_y, phi_x, matches("^phi_(.*)_power$"), matches("^phi1_(.*)_power$"))) %>%
  pivot_longer(phi_adl_power:phi1_adl33_power, names_to = "model", values_to = "phi1_power") %>%
  mutate(model = gsub("^ph.*_(.*)_power", "\\1", model))
phi2_bias <- data %>%
  select(c(time, phi_y, phi_x, matches("^bias_(.*)_phi2$"))) %>%
  pivot_longer(bias_adl22_phi2:bias_adl33_phi2, names_to = "model", values_to = "phi2_bias") %>%
  mutate(model = gsub("^bias_(.*)_phi2", "\\1", model))
phi2_type1 <- data %>%
  select(c(time, phi_y, phi_x, matches("^phi2_(.*)_t1$"))) %>%
  pivot_longer(phi2_adl22_t1:phi2_adl33_t1, names_to = "model", values_to = "phi2_type1") %>%
  mutate(model = gsub("^ph.*_(.*)_t1", "\\1", model))
phi3_bias <- data %>%
  select(c(time, phi_y, phi_x, matches("^bias_(.*)_phi3$"))) %>%
  pivot_longer(bias_adl33_phi3, names_to = "model", values_to = "phi3_bias") %>%
  mutate(model = gsub("^bias_(.*)_phi3", "\\1", model))
phi3_type1 <- data %>%
  select(c(time, phi_y, phi_x, matches("^phi3_(.*)_t1$"))) %>%
  pivot_longer(phi3_adl33_t1, names_to = "model", values_to = "phi3_type1") %>%
  mutate(model = gsub("^ph.*_(.*)_t1", "\\1", model))
b1_bias <- data %>%
  select(c(time, phi_y, phi_x, matches("^bias_(.*)_b1$"))) %>%
  pivot_longer(bias_adl_b1:bias_adl33_b1, names_to = "model", values_to = "b1_bias") %>%
  mutate(model = gsub("^bias_(.*)_b1", "\\1", model))
b1_type1 <- data %>%
  select(c(time, phi_y, phi_x, matches("^b1_(.*)_t1$"))) %>%
  pivot_longer(b1_adl_t1:b1_adl33_t1, names_to = "model", values_to = "b1_type1") %>%
  mutate(model = gsub("^b1_(.*)_t1", "\\1", model))
b1_power <- data %>%
  select(c(time, phi_y, phi_x, matches("^b1_(.*)_power$"))) %>%
  pivot_longer(b1_adl_power:b1_adl33_power, names_to = "model", values_to = "b1_power") %>%
  mutate(model = gsub("^b1_(.*)_power", "\\1", model))
b2_bias <- data %>%
  select(c(time, phi_y, phi_x, matches("^bias_(.*)_b2$"))) %>%
  pivot_longer(bias_adl_b2:bias_adl33_b2, names_to = "model", values_to = "b2_bias") %>%
    mutate(model = gsub("^bias_(.*)_b2", "\\1", model))
b2_type1 <- data %>%
  select(c(time,phi_y, phi_x, matches("^b2_(.*)_t1$"))) %>%
  pivot_longer(b2_adl_t1:b2_adl33_t1, names_to = "model", values_to = "b2_type1") %>%
  mutate(model = gsub("^b2_(.*)_t1", "\\1", model))
b2_power <- data %>%
  select(c(time, phi_y, phi_x, matches("^b2_(.*)_power$"))) %>%
  pivot_longer(b2_adl_power:b2_adl33_power, names_to = "model", values_to = "b2_power") %>%
  mutate(model = gsub("^b2_(.*)_power", "\\1", model))
b3_bias <- data %>%
  select(c(time, phi_y, phi_x, matches("^bias_(.*)_b3$"))) %>%
  pivot_longer(bias_adl22_b3:bias_adl33_b3, names_to = "model", values_to = "b3_bias") %>%
  mutate(model = gsub("^bias_(.*)_b3", "\\1", model))
b3_type1 <- data %>%
  select(c(time, phi_y, phi_x, matches("^b3_(.*)_t1$"))) %>%
  pivot_longer(b3_adl22_t1:b3_adl33_t1, names_to = "model", values_to = "b3_type1") %>%
  mutate(model = gsub("^b3_(.*)_t1", "\\1", model))
b4_bias <- data %>%
  select(c(time, phi_y, phi_x, matches("^bias_(.*)_b4$"))) %>%
  pivot_longer(bias_adl33_b4, names_to = "model", values_to = "b4_bias") %>%
  mutate(model = gsub("^bias_(.*)_b4", "\\1", model))
b4_type1 <- data %>%
  select(c(time, phi_y, phi_x, matches("^b4_(.*)_t1$"))) %>%
  pivot_longer(b4_adl33_t1, names_to = "model", values_to = "b4_type1") %>%
  mutate(model = gsub("^b4_(.*)_t1", "\\1", model))
lrm_bias <- data %>%
  select(c(time, phi_y, phi_x, matches("^bias_(.*)_lrm$"))) %>%
  pivot_longer(bias_adl_lrm:bias_adl33_lrm, names_to = "model", values_to = "lrm_bias") %>%
  mutate(model = gsub("^bias_(.*)_lrm$", "\\1", model))
lrm_type1 <- data %>%
  select(c(time, phi_y, phi_x, matches("^lrm_(.*)_t1$"))) %>%
  pivot_longer(lrm_adl_t1:lrm_adl33_t1, names_to = "model", values_to = "lrm_type1") %>%
  mutate(model = gsub("^lrm_(.*)_t1", "\\1", model))
lrm_power <- data %>%
  select(c(time, phi_y, phi_x, matches("^lrm_(.*)_power$"))) %>%
  pivot_longer(lrm_adl_power:lrm_adl33_power, names_to = "model", values_to = "lrm_power") %>%
  mutate(model = gsub("^lrm_(.*)_power", "\\1", model))
ftest_y_t1 <- data %>%
  dplyr::select(c(time,phi_y, phi_x, matches("^ftest_y_(.*)_t1$"))) %>%
  pivot_longer(ftest_y_adl33_t1:ftest_y_adl22_t1, names_to = "model", values_to = "ftest_y_t1") %>%
  mutate(model = gsub("^ftest.*_(.*)_t1", "\\1", model))
ftest_x_t1 <- data %>%
  dplyr::select(c(time,phi_y, phi_x, matches("^ftest_x_(.*)_t1$"))) %>%
  pivot_longer(ftest_x_adl33_t1:ftest_x_adl22_t1, names_to = "model", values_to = "ftest_x_t1") %>%
  mutate(model = gsub("^ftest.*_(.*)_t1", "\\1", model))
ftest_y_power <- data %>%
  dplyr::select(c(time,phi_y, phi_x, matches("^ftest_y_(.*)_power$"))) %>%
  pivot_longer(ftest_y_adl33_power:ftest_y_adl22_power, names_to = "model", values_to = "ftest_y_power") %>%
  mutate(model = gsub("^ftest.*_(.*)_power", "\\1", model))
ftest_x_power <- data %>%
  dplyr::select(c(time,phi_y, phi_x, matches("^ftest_x_(.*)_power$"))) %>%
  pivot_longer(ftest_x_adl33_power:ftest_x_adl22_power, names_to = "model", values_to = "ftest_x_power") %>%
  mutate(model = gsub("^ftest.*_(.*)_power", "\\1", model))

#----Ly Bias----
phi1_bias$model <- factor(phi1_bias$model, levels = c("gecm", "adl", "adl22", "adl33"), labels = c("gecm", "adl", "adl22", "adl33"))
ggplot(phi1_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi1_bias), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Bias", x = "T", shape = "", title = expression("Bias in the coefficient of y"[t-1])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +  
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_phi1_bias.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----Ly Type 1----
phi1_type1$model <- factor(phi1_type1$model, levels = c("gecm", "adl", "adl22", "adl33"), labels = c("gecm", "adl", "adl22", "adl33"))
ggplot(phi1_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi1_type1), position = position_dodge(width = 0.65)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Type 1 error rate in the coefficient of y"[t-1])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_phi1_type1.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----Ly Coverage----
phi1_type1 %>%
  mutate(phi1_coverage = 1 - phi1_type1) %>%
  ggplot(aes(shape = model)) +
    geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
    geom_point(aes(x = time, y = phi1_coverage), position = position_dodge(width = 0.65)) +
    theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
    labs(y = "Coverage Probability", x = "T", shape = "", title = expression("Coverage probability of the coefficient of y"[t-1])) + ylim(0, 1) +
    theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
          panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
          panel.grid.minor.y = element_blank()) + 
    scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
    scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                       values = c(3, 4, 5, 6),
                       labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
    annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
    annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
    annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
    annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
    annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
    annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
    annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
    facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_phi1_coverage.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----Ly Power----
phi1_power$model <- factor(phi1_power$model, levels = c("gecm", "adl", "adl22", "adl33"), labels = c("gecm", "adl", "adl22", "adl33"))
ggplot(phi1_power, aes(shape = model)) +
  geom_hline(yintercept = 1, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi1_power), position = position_dodge(width = 0.65)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Power", x = "T", shape = "", title = expression("Power in the coefficient of y"[t-1])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_phi1_power.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----L2y Bias----
phi2_bias$model <- factor(phi2_bias$model, levels = c("adl22", "adl33"), labels = c("adl22", "adl33"))
ggplot(phi2_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi2_bias), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Bias", x = "T", shape = "", title = expression("Bias in the coefficient of y"[t-2])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)"))  +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_phi2_bias.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----L2y Type 1----
phi2_type1$model <- factor(phi2_type1$model, levels = c("adl22", "adl33"), labels = c("adl22", "adl33"))
ggplot(phi2_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi2_type1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Type 1 error rate in the coefficient of y"[t-2])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +  
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)"))  +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_phi2_type1.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----L2y Coverage----
phi2_type1 %>%
  mutate(phi2_coverage = 1 - phi2_type1) %>%
ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi2_coverage), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x = "T", shape = "", title = expression("Coverage probability of the coefficient of y"[t-2])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +  
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)"))  +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_phi2_coverage.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----L3y Bias----
ggplot(phi3_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi3_bias), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Bias", x = "T", shape = "", title = expression("Bias in the coefficient of y"[t-3])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) +
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl33"),
                     values = c(6),
                     labels = c("ADL(3,3)")) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x))) 
ggsave("ADL_phi3_bias.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----L3y Type 1----
ggplot(phi3_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi3_type1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Type 1 error rate in the coefficient of y"[t-3])) + ylim(0, 1) + 
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) +
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl33"),
                     values = c(6),
                     labels = c("ADL(3,3)")) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x))) 
ggsave("ADL_phi3_type1.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----L3y Coverage----
phi3_type1 %>%
  mutate(phi3_coverage = 1 -phi3_type1) %>%
ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi3_coverage), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x = "T", shape = "", title = expression("Coverage probability of the coefficient of y"[t-3])) + ylim(0, 1) + 
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) +
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl33"),
                     values = c(6),
                     labels = c("ADL(3,3)")) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x))) 
ggsave("ADL_phi3_coverage.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----x Bias----
b1_bias$model <- factor(b1_bias$model, levels = c("gecm", "adl", "adl22", "adl33"), labels =  c("gecm", "adl", "adl22", "adl33"))
ggplot(b1_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b1_bias), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Bias", x = "T", shape = "", title = expression("Bias in the coefficient of x"[t])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x))) 
ggsave("ADL_b1_bias.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----x Type 1----
b1_type1$model <- factor(b1_type1$model, levels = c("gecm", "adl", "adl22", "adl33"), labels =  c("gecm", "adl", "adl22", "adl33"))
ggplot(b1_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b1_type1), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Type 1 error rate in the coefficient of x"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x))) 
ggsave("ADL_b1_type1.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----x Coverage----
b1_type1 %>%
  mutate(b1_coverage = 1 - b1_type1) %>%
ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b1_coverage), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x = "T", shape = "", title = expression("Coverage probability of the coefficient of x"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x))) 
ggsave("ADL_b1_coverage.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----x Power----
b1_power$model <- factor(b1_power$model, levels = c("gecm", "adl", "adl22", "adl33"), labels =  c("gecm", "adl", "adl22", "adl33"))
ggplot(b1_power, aes(shape = model)) +
  geom_hline(yintercept = 1, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b1_power), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Power", x = "T", shape = "", title = expression("Power in the coefficient of x"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x))) 
ggsave("ADL_b1_power.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----Lx Bias----
b2_bias$model <- factor(b2_bias$model, levels = c("gecm", "adl", "adl22", "adl33"), labels =  c("gecm", "adl", "adl22", "adl33"))
ggplot(b2_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b2_bias), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Bias", x = "T", shape = "", title = expression("Bias in the coefficient of x"[t-1])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x))) 
ggsave("ADL_b2_bias.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----Lx Type 1----
b2_type1$model <- factor(b2_type1$model, levels = c("gecm", "adl", "adl22", "adl33"), labels =  c("gecm", "adl", "adl22", "adl33"))
ggplot(b2_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b2_type1), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Type 1 error rate in the coefficient of x"[t-1])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)"))  +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_b2_type1.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----Lx Coverage----
b2_type1 %>%
  mutate(b2_coverage = 1 - b2_type1) %>%
ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b2_coverage), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x = "T", shape = "", title = expression("Coverage probability of the coefficient of x"[t-1])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)"))  +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_b2_coverage.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----Lx Power----
b2_power$model <- factor(b2_power$model, levels = c("gecm", "adl", "adl22", "adl33"), labels =  c("gecm", "adl", "adl22", "adl33"))
ggplot(b2_power, aes(shape = model)) +
  geom_hline(yintercept = 1, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b2_power), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Power", x = "T", shape = "", title = expression("Power in the coefficient of x"[t-1])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) + 
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)"))  +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_b2_power.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----L2x Bias----
ggplot(b3_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b3_bias), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Bias", x = "T", shape = "", title = expression("Bias in the coefficient of x"[t-2])) + ylim(-1, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)"))  +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_b3_bias.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----L2x Type 1----
ggplot(b3_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b3_type1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Type 1 error rate in the coefficient of x"[t-2])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)"))  +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_b3_type1.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----L2x Coverage----
b3_type1 %>%
  mutate(b3_coverage = 1 - b3_type1) %>%
ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b3_coverage), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x = "T", shape = "", title = expression("Coverage probability of the coefficient of x"[t-2])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)"))  +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_b3_coverage.pdf", width = 25, height = 10, units = "in", dpi = 300)


#----L3x Bias----
ggplot(b4_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b4_bias), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Bias", x = "T", shape = "", title = expression("Bias in the coefficient of x"[t-3])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl33"),
                     values = c(6),
                     labels = c("ADL(3,3)"))  +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_b4_bias.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----L3x Type 1----
ggplot(b4_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b4_type1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Type 1 error rate in the coefficient of x"[t-3])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl33"),
                     values = c(6),
                     labels = c("ADL(3,3)"))  +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_b4_type1.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----L3x Coverage----
b4_type1 %>%
  mutate(b4_coverage = 1 - b4_type1) %>%
ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b4_coverage), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x = "T", shape = "", title = expression("Coverage probability of the coefficient of x"[t-3])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl33"),
                     values = c(6),
                     labels = c("ADL(3,3)"))  +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_b4_coverage.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----LRM Bias----
lrm_bias$model <- factor(lrm_bias$model, levels = c("gecm", "adl", "adl22", "adl33"), labels = c("gecm", "adl", "adl22", "adl33"))
ggplot(lrm_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = lrm_bias), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Bias", x = "T", shape = "", title = expression("Bias in the LRM of x"[t])) + ylim(-20, 20) + 
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c( 3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)"))  +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_lrm_bias.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----LRM Type 1----
lrm_type1$model <- factor(lrm_type1$model, levels = c("gecm", "adl", "adl22", "adl33"), labels = c("gecm", "adl", "adl22", "adl33"))
ggplot(lrm_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = lrm_type1), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Type 1 error rate in the LRM of x"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)"))  +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_lrm_type1.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----LRM Coverage----
lrm_type1 %>%
  mutate(lrm_coverage = 1 - lrm_type1) %>%
ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = lrm_coverage), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x = "T", shape = "", title = expression("Coverage probability of the LRM of x"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)"))  +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_lrm_coverage.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----LRM Power----
lrm_power$model <- factor(lrm_power$model, levels = c("gecm", "adl", "adl22", "adl33"), labels = c("gecm", "adl", "adl22", "adl33"))
ggplot(lrm_power, aes(shape = model)) +
  geom_hline(yintercept = 1, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = lrm_power), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Power", x = "T", shape = "", title = expression("Power in the LRM of x"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("gecm", "adl", "adl22", "adl33"),
                     values = c(3, 4, 5, 6),
                     labels = c("GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)"))  +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_lrm_power.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----F-test, Lags of y Type 1 (sum = true phi_1)----
ftest_y_t1 %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_y_t1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x= "T", shape = "", title = expression("Type 1 error rate for a F-test on the coefficients of the lags of y"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_ftest_y_t1.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----F-test, Lags of y Coverage (sum = true phi_1)----
ftest_y_t1 %>%
  mutate(ftest_y_coverage = 1 - ftest_y_t1) %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_y_coverage), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x= "T", shape = "", title = expression("Coverage probability for a F-test on the coefficients of the lags of y"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_ftest_y_coverage.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----F-test, Lags of y Power (sum = 0 for non-0 phi_1)----
ftest_y_power %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 1, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_y_power), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Power", x= "T", shape = "", title = expression("Power for a F-test on the coefficients of the lags of y"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_ftest_y_power.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----F-test, Lags of x Type 1 (sum = tru beta + gamma)----
ftest_x_t1 %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_x_t1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x= "T", shape = "", title = expression("Type 1 error rate for a F-test on the coefficients of x"[t]~"and its lags")) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_ftest_x_type1.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----F-test, Lags of x Coverage (sum = tru beta + gamma)----
ftest_x_t1 %>%
  mutate(ftest_x_coverage = 1 - ftest_x_t1) %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_x_coverage), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x= "T", shape = "", title = expression("Coverage probability for a F-test on the coefficients of x"[t]~"and its lags")) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_ftest_x_coverage.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----F-test, Lags of x Power (sum = 0)----
ftest_x_power %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 1, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_x_power), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Power", x= "T", shape = "", title = expression("Power for a F-test on the coefficients of x"[t]~"and its lags")) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~phi_x, labeller = label_bquote(phi[y] == .(phi_y), phi[x] == .(phi_x)))
ggsave("ADL_ftest_x_power.pdf", width = 25, height = 10, units = "in", dpi = 300)
